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ABSTRACT 

We continue exploration of the Boltzmann scheme started in Banerjee and 
Ghosh (2007, henceforth Paper I) for studying the evolution of compact-binary 
populations of globular clusters, introducing in this paper our method of handling 
the stochasticity inherent in dynamical processes of binary formation, destruc- 
tion and hardening in globular clusters. We describe these stochastic processes 
as Wiener processes, whereupon the Boltzmann equation becomes a stochastic 
partial differential equation, the solution of which requires the use of ltd calcu- 
lus (this use being the first, to our knowledge, in this subject), in addition to 
ordinary calculus. We focus on the evolution of (a) the number of X-ray binaries 
NxB in globular clusters, and (b) the orbital-period distribution of these binaries. 
We show that, although the details of the fluctuations in the above quantitities 
differ from one "realization" to another of the stochastic processes, the general 
trends follow those found in the continuous-limit study of Paper I, and the aver- 
age result over many such realizations is close to the continuous-limit result. We 
investigate the dependence of Nxb found by these calculations on two essential 
globular-cluster parameters, namely, the star-star and star-binary encounter-rate 
parameters P and 7, for which we had coined the name Verbunt parameters in 
Paper I. We compare our computed results with those from CHANDRA obser- 
vations of Galactic globular clusters, showing that the expected scalings of Nxb 
with the Verbunt parameters are in good agreement with the observed ones. We 
indicate what additional features can be incorporated into the scheme in future, 
and how more elaborate problems can be tackled. 



Subject headings: globular clusters: general — binaries: close — X-rays: binaries 
— methods: numerical — stellar dynamics — scattering 
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Introduction 



In this series of papers, we are studying the evolution of compact-bina ry populations of 
globul ar clusters with the aid of a Boltzmann scheme which we introduced in lBanerjee fc Ghosh 
( 120071 ). henceforth Paper I. This scheme follows compact-binary evolution as a result of both 
(a) those processes which determine compact-binary evolution in isolation {i.e., outside glob- 
ular clusters), e.g., angular momentum loss by gravitational radiation and magnetic braking, 
as also orbital evolution due to mass transfer, and (b) those processes which arise from en- 
counters of compact binaries with t he dense ste l lar background in globul ar clusters, e.g., 
collisional hardening (IHeggie I Il975l : IShull 1 1 19791 : iBanerjee fc Ghosh 1 12006| ) , binary forma- 



tion t hrough tidal capture and exchange processes, and binary destruction (IFabian et.al. 



1975 



Press fc Teukolskv Ill977| : Ihee fc Ostriker I ll986l: loi Stefano fc Rappaport Ill992l . 



1994 



Spitzer 1 119871 : iHut &: Bahcall I Il983l ). We treat all of the above processes simultaneously 
through our Boltzmann scheme, the aim being to see their combined effect on the compact- 
binary population as a whole, in particular on the evolution of (a) the total number of X-ray 
binaries as the formation and destruction processes continue to operate, and, (b) the orbital- 
period distribution of the population. As stressed in Paper I, our scheme is the original 
Boltzmann one (not the Fokker-Planck reduction of it), which, by definition, is capable of 
handling both the combined small effects of a large number of frequent, weak, distant encoun- 
ters and the individual large effects of a small number of rare, strong, close encounters on 
the same footing. We note here that, although Monte Carlo Fokker-Planck approaches were 
normally thought to be capable of handling only th e former effects, sch e mes for including 



the l atter have been proposed and studied recently (IFregeau et.al I l2003l : iFregeau fc Rasio 
2007h . 



In Paper I, we studied the problem in the continuous limit, wherein we used continu- 
ous representations for both kinds of processes described above, i.e., those of category (a) 
above, which are inherently continuous, and also those of category (b), which are inherently 
stochastic. For the latter category, therefore, we used the continuous limit of the above 
stochastic processes, wherein the probability or cross-section of a particular such process 
happening with a given set of input and output variables was treated as a continuous func- 
tion of these variables. These cross- sect ions were, of course, those that had been determined 
from ex tensive numerical experinients with two-body and three-body encounters performed 
earlier Jfleggie. Hut fc McMillan] Il996l : [Porteries Zwart et.al lE997bf ). 



In this paper, we address the next question, namely, how is the inherent stochasticity of 
the processes of category (b) to be introduced into our scheme, to be handled simultaneously 
with the inherently continuous nature of those of category (a)? As stressed in Paper I, this 
step is of great importance, since it is a simultaneous operation of the above continuous and 
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stochastic processes in globular clusters that leads to the observed properties of compact- 
binary populations in them. To this end, we introduce stochasticity into our Boltzmann 
study in this paper in the following way. For a first look, we consider the rates of the 
processes of category (b) as randomly fluctuating about the mean rates described in Paper 
I, while those of the processes of category (a) remain continuous, as before. We model 
these fluctuations as a Wiener process (see Appendix A and references therein), which is the 
mathematical description of Brownian motion. 

With this prescription, the Boltzmann equation governing the evolution of the distribu- 
tion function n{a,t) of compact binaries in time t and orbital radius a becomes a stochastic 
partial differential equation (henceforth SPDE), instead of the ordinary partial differential 
equation (henceforth OPDE) which it was in the continuous limit. We handle the solu- 
tion of this SPDE wi t h the aid of t echniques develo ped largely during the last fifteen years 



(IKloeden et.al 1 1 19941 : I Gains 1 1 19951 : l0ksendal 1 120041 ) . These techniques involve the use of 



the ltd calculus (see Appendix B and references therein), instead of ordinary calculus, for 
handling the stochastic terms. 

Our results show that the full solutions with stochasticity included have fluctuations 
which vary from one "realization" to another of the stochastic processes, as expected. How- 
ever, the full results show trends which generally follow those in the continuous limit. Fur- 
thermore, the average result over many realizations comes very close to the continuous limit, 
showing the importance of the latter limit for understanding mean trends. On the other 
hand, understanding fluctuations in a typical full run is also very important, as this gives 
us a first idea of the magnitude of fluctuations we can expect in the data on X-ray binaries 
in globular clusters as a result of the stochastic processes, as also the expected trends in 
the fluctuations with the essential globular-cluster parameters, e.g., the Verbunt parameters 
introduced in Paper I (also see below). 

Comparison of our computed trends in the number Nxb of X-ray binaries in Galactic 
globular clusters with the Verbunt parameters on the one hand, with observed trends in 
recent CHANDRA data on Galactic globular clusters on the other, shows that our full results 
are in good agreement with observation. We have thus constructed a straightforward, very 
inexpensive scheme for following the evolution of compact-binary populations in globular 
clusters, including essential, fluctuating, encounter processes that are thought to operate in 
such clusters, as also those continuous processes which operate in isolated binaries and so 
apply here as well. We can also follow the evolution of Nxb, as also that of the orbital-period 
distribution of compact binaries in globular clusters. For the latter study, however, proper 
modeling of stellar-evolutionary effects still remains to be done for parts of the parameter 
space, as explained in Paper I, and as discussed in Sec. [51 
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In Sec. |2l we briefly review the continuous-limit results of Paper I, in order to put 
the results of this paper in their proper context. We give only the essentials here, citing 
figures in Paper I for detailed results. In Sec. [3|, we introduce stochasticity explicitly through 
our prescription, explaining the details of Wiener processes and the Ito calculus in the 
Appendices. We describe our generalization of the Lax-Wendorff scheme, introduced in 
Paper I, to handle the solution of the SPDE which the Boltzmann equation has become 
now. In Sec. HI we describe the results of our full calculations including stochasticity, and 
compare these with the continuous-limit results of Paper I. In Sec. 14.31 we compare our full 
results with observations. Finally, In Sec. [5|, we discuss our results, putting them in the 
context of previous studies in the subject, and indicating some additional physical effects to 
be included by stages in future versions of our scheme, as well as some future problems to 
be tackled. 



2. Brief Review of Continuous Limit 

In order to put the stochastic studies of this paper in their proper context, we review in 
this section the essentials of the continuous- limit studies of Paper I which form this context. 
In the latter limit, the Boltzmann description works in terms of appropriate mean values of 
the variables and parameters which are actually stochastic. Accordingly, the above compact- 
binary distribution function n{a, t) is replaced by its mean value n{a, t), and the Boltzmann 
equation has the form: 

^=Ma)~n(a,ma)-^7(a). (1) 

Here, R{a) is the mean formation rate (per unit binary radius) of compact binaries of radius 
a, D{a) is the mean destruction rate per binary and /(a) is the mean shrinkage rate d of a 
compact binary of radius a, as described in Paper I. 



2.1. Mean rates 

The mean shrinkage or "hardening" rate /(a) = a has been given in Fig. 2 of Paper I as 
a function of a, describing the situation as the compact binary goes from its widest, pre-X- 
ray-binary (PXB) phase to Roche lobe contact, and continues through the mass-transferring 
X-ray-binary (XB) phase. As shown there, collisional hardening, i.e., that due to encounters 
between the binary and the stellar background of the globular cluster, dominates at large a, 
while hardening by gravitational radiation and magnetic braking dominates at small a. The 
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relative orbit shrinkage rate a /a scales roughly as a at large orbital radii, passes through a 
minimum at a critical separation where the gravitational radiation shrinkage rate, scaling 
as d/a ~ a~^, takes over. Magnetic braking also contributes at small radii, but Roche lobe 
contact also occurs at roughly the same point, whereupon the angular-momentum transfer 
associated with mass transfer in the XB phase dominates the orbit-change rate, and a has 
a very weak dependence on a during this phase. Detailed quantitative expressions for the 
above rates are given in Paper I, to which we refer the reader, recording here only that the 
orbit-shrinkage rate is given in terms of the angular-momentum loss rate by 

a J fric rhx , ^ 

- = 2--2— -2— , 2 
a J rric mx 

where mx is mass of the degenerate star in the compact binary, and rric is that of its (low- 
mass) companion, and the total angular-momentum loss rate can be written in terms of its 
components as: 

J (a) = j= jGw{a) + ]MB{a) +jcoii{a)- (3) 

Here, the subscripts 'GW, 'MB' and 'coll' respectively stand for gravitational radiation, 
magnetic braking, and collisional hardening. 

The other essential mean rates are those of compact binary formation and destruction, 
R{a) and D{a), respectively. Consider first the formation of compact binaries with degen- 
erate primaries (white dwarfs or neutron stars) and low-mass companions, such as we are 
interested in this series of papers, in globular cluster (henceforth GC) cores. The two relevant 
dynamical processes are (i) tidal capture (tc) of a degenerate, compact star by an ordinary, 
low-mass star, and (ii) an exchange encounter (exl) between such a compact star and a 
binary of two ordinary low-mass stars in the GC, wherein the compact star replaces one of 
the binary members. Accordingly, the total mean rate of formation of compact binaries per 
unit binary radius, R[a), consists of the above mean tc rate rtc{a) and mean exl rate Vexiia)'- 

R{a) =rtcia) + r,^i{a). (4) 

The above mean rates are shown as functions of a in Fig. 3 of Paper I, and detailed mathe- 
matical expressions for them are also given in that reference, which we do not repeat here. 
The mean tidal-capture rate is nearly constant for a < 5Rq, and decreases rapidly at larger 
a. Extensive discussion of various issues related to tidal capture and of the current status 
of our understanding of this process are also given in Paper I, to which we refer the reader. 
Further discussion of previous studies of tidal capture in this problem are given in Sec. O 
The mean exchange (exl) rate is roughly constant over the range of radii of interest here, for 
the widely-adopted radius distribution of primordial binaries, viz., a uniform distribution in 
In a, which we adopt throughout our work. 
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Consider now the destruction of compact binaries with degenerate primaries and low- 
mass companions, which can occur in the following two ways. First, an encounter with a 
star which has a relative speed higher than an appropriate critical speed can lead to its 
dissociation (dss). Second, in an exchange encounter (ex2) of this binary with a compact 
star, the latter can replace the low-mass companion in the binary, forming a double compact- 
star binary consisting of two neutron stars, two white dwarfs, or a neutron star and a white 
dwarf (all with masses mx ~ IAMq in our model: see Paper I). This destroys the binary 
as an X-ray source (as accretion is not possible in such a system), and so takes it out of 
reckoning in our studjfl. The total mean destruction rate D{a) per binary is thus the sum 
of the above mean dss and ex2 rates: 



The above mean rates are shown as functions of a in Fig. 3 of Paper I, and detailed math- 
ematical expressions for them are given in that paper, which, once again, we do not repeat 
here. The mean dissociation rate is negligible below a critical radius Oc corresponding to 
the above critical speed, rises extremely sharply above Oc at first, and eventually scales as 
for a ^ ttc- By contrast, the mean exchange (ex2) rate is roughly oc a and dominates 
the destruction processes completely at all orbital radii relevant to our study, dissociation 
becoming important only for very soft binaries (a > lOOOi?©, say), which are of little interest 
to us here. 



We now summarize the essentials of our continuous-limit results in Paper I. We com- 
puted evolution of the compact-binary distribution in this limit, showing the surface n{a,t) 
explicitly in three dimensions in Fig. 5 of Paper I, corresponding to representative GC pa- 
rameters (rather similar to those of the well-known Galactic cluster 47 Tuc). The surface 
evolved smoothly, with the compact-binary population growing predominantly at shorter 
radii (a < IORq, say). We showed that, starting with a small number of binaries at t = 
following various distributions, we obtained at times ~ Gyr or longer a distribution which 
was independent of the initial conditions, determined entirely by the dynamical processes 
of formation and destruction, and by the various hardening processes summarized above. 
We clarified the nature of the distribution and its evolution by displaying slices through the 



^Note that it is essentially impossible for one of the compact stars in such a double-compact system to 
be re-exchanged with an ordinary star in a subsequent exchange encounter, since the average mass of a 
background GC star (taken as m/ = O.6A/0 in our work) is much less than the above value of mx- 




(5) 



2.2. 



Results 
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above surface at various points along time axis and a-axis, shown in Figs. 6 and 7 of Paper 
I. The former figure showed that the profile n{a) increased with time, roughly preserving 
its profile for t > 1.5 Gyr, this profile consisting of a roughly uniform distribution at small 
orbital radii, a < 6Rq, say, and a sharp fall-off at larger radii. The latter figure showed that 
n{a) at a given a increased with time and approached saturation on a timescale 6 — 12 Gyr, 
the timescale being longer at smaller values of a. 

For comparison with crucial X-ray observations of Galactic GCs, we computed the total 
number of XBs Nxb in a GC at any time, obtained by integrating n{a,t) over the range of 
a relevant for XBs, viz., apm ^ a (11, where apm is the value of a corresponding to the 
period minimum P ^ 80 minutes, and is the value of a at the first Roche lobe contact 
and onset of mass transfer, as explained above: 



NxB{t) 



n{a, t)da. 



(6) 



Taking an evolutionary time ~ 8 Gyr as representative, we determined Nxb at this point 
in time, and studied its dependence on the Verbunt parameters F and 7 that describe the 
essential dynamical properties of globular clusters in this context, as explained in Paper I, 
showing our results as the computed Nxb^Xi 7) surface in three dimensions in Fig. 8 of that 
paper. 

The Verbunt parameters F and 7 hav e been introduced in Paper I. Following pione ering 
suggestions by Verbunt and co-authors (I Verbunt &: Hut Ml 98 71 : IVerbunt et.al I Il989l ). the 
crucial impo rtance of these param eters in GC dynamics has been lucidly summarized recently 
by Verbunt ( Verbunt II20031 . 120061 ). and the importance of the parameter F for scaling between 
different G Cs has been emphasized in a pion eering study of the production of recycled pulsars 
in GCs by iDi Stefano fc Rappaport I (119921 ). Briefly, the first parameter F is the two-body 
stellar encounter rate, which scales with p^r^/vc, and occurs naturally in the rates of all 
two-body processes, where the standard GC core variables are the average stellar density p, 
the velocity dispersion Vc, and the core radius Tc- In fact, we defined F as 



r ^ ^ oc p'/^l 



(7) 



in Paper I. Note that the last scaling in the above equation holds only for virialized cores, 
where the scaling Vc oc p^l'^r^. can be applied. The second parameter is a measure of the rate 
of encounter between binaries and single stars in the GC, the rate normally used being the 
encounter rate 7 of a single binary with the stellar background, with the understanding that 
the total rate of binary-single star encounter in the cluster will be oc 717. We defined 7 in 
Paper I as 

7 = ^ oc p'l\z\ (8) 
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where the last scahng holds, again, only for virialized cores. 

We compared the results of our above computations with the systema tics of recent 



observations of X-ray binaries in Galactic globular clusters (iPooley et.al II2003I ). as displayed 
in Figs. 9 and 10 of Paper I. We showed that the computed total number Nxb of XBs 
expected in a globular cluster scaled in a characteristic way with the Verbunt parameters. 
The qu alitative nature of this sc aling was rather similar to that found in our earlier "toy" 



model (jBanerjee fc Ghosh II2006I ). although details were different. Nxb scaled with F (which 
is a measure of the dynamical formation rate of compact binaries, as above) and, at a 
given F, Nxb decreased with increasing 7 (which is a measure of the rate of destruction of 
these binaries by dynamical processes) at large values of 7 , as shown in Fig. 9 of Paper 
I. This rough scaling could be expressed as Nxb oc Tgi^j), where the "universal" function 
(7(7) of 7 (except for a spurious feature at low values of 7 which we explained in Paper I) 
decreased monotonically with increasing 7, reflecting the increasing strength of dynamical 
binary-destruction processes with increasing 7. 

We further demonstrated that these computed trends with the Verbunt parameters 
compared very well with the observed trends in above X-ray data by showing in Fig. 10 of 
Paper I the contours of constant Nxb in the F — 7 plane, and overplotting on the positions 
of the Galactic GCs with observed X-ray binaries in them. We showed that the trend in the 
observed Nxb values generally followed the contours, with one exception. This provided us 
with a first indication of the basic ways in which dynamical binary formation and destruction 
processes work in GCs, and encouraged us to build more "realistic" models by introducing 
stochastic effects explicitly, to which this paper is devoted. 



3. Introducing Stochasticity 

In order to study the behavior of the inherently stochastic terms in the full Boltzmann 
equation 

= Ria, t) - nia, t)D{a, t) - ^^fia, t), (9) 

we must explicitly include stochastic, fluctuating parts in these terms, in addition to their 
mean values studied in Paper I, as above. We do so by expressing the above rates R{a,t), 
D{a,t), and f{a,t) as their earlier mean values R{a), D{a) and /(a), augmented by fluctu- 
ating components as below: 

i?(a, t) = Ria) + Cat. + 
D(a,t) = D(a) + C.2 + CLs. ) (10) 



/(a,t) = /(a) + C 



t 

a coll 
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Here, is the random fluctuation rate of events of type X from tlieir mean rates, and X 
= tc, exl, ex2, dss, coll by turn, these notations having been introduced above. In general, 
(^j^ is a function of both a and t, of course. 

The crucial question is that of modeling appropriately. In this introductory work, 
we use the standard normally-distributed model 

= Sx{a)v\ (11) 

where Sx{a) is the variance of C*^^ at a given a and rj^s at each t are independent random 
numbers distributed in a standard normal distribution. This separable form is appropriate 
since the dynamical processes of binary formation and destruction at a given value of a are 
inherently independent of those at other values of a. The "flow" or "current" of binaries 
from larger to smaller values of a due to the hardening described above and in Paper I does 
not affect this independence, but merely changes the number of binaries in an infinitesimal 
interval of a around a given value of a at a given instant t, which is automatically taken into 
account by the Boltzmann equation (also see below). Indeed, the hardening process itself 
has this independence, viz., that its rate at a given value of a is independent of that at other 
values of a, and so is separable in the same way. By contrast, the number distribution n{a, t) 
of the binaries cannot be written in this form, since, at a particular a, it is determined both 
by the binary formation and destruction rates at that a, and by the rates of binary arrival 
from (and also departure to) other values of a due to hardening, as described above. All of 
this is, of course, automatically included in the Boltzmann equation by definition. 

The essence of the physics of these fluctuations is contained in the adopted model for r/*. 
By adopting a normally-distributed variation, we are, in effect, considering a Wiener process 
(see Appendix A and references therein), which is the standard mathematical description of 
Brownian motion. In other words, we are studying a situation wherein the variations in the 
above dynamical rates about their respective mean values constitute a Brownian motion. 
We return to Wiener processes later in more detail. 



3.1. Variances of stochastic-process rates 



How do we estimate the variance of a stochastic process of type X whose mean value 
is Rx{cl)1 To answer this question, consider first how it is addressed in Monte Carlo 
simulations, which have bee n performed in this subject b y several authors ( se e, e.g., 



Sigurdsson &: Phinney I (119931 ) . iPortegies Zwart et.al I (Il997al ). or lFregeau et.al I (120031 )). These 



works have uniformly used the so-called rejection method for determining whether an event 
of a given type occurs in a given time interval or not. The method works as follows. 
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For events of type X, if the mean rate of event occurrence is Rx-, then the timescale for 
occurrence of such events is 

Atx = =- (12) 

Hence, during a time step At < Atx, the quantity px = Rx^t < 1 is the expected mean 
number of events during this interval, px < 1 can also be interpreted as t he probability of 



occurrence of an event X within this time step (jPortegies Zwart et.al Ill997al ). and the actual 



number of such events within At will then follow a binomial distribution with the following 
mean and variance: 



mean = Rx{a)At 
variance = Sj^{a)At'^ = Rxia)At{l - Rxia)At). 



(13) 



Note that the above variance depends on a, since the mean rates depend on a. When several 
different types of events are considered simultaneously, as in the present problem, we must, of 
course, so choose At that it is shorter than the shortest event-occurrence timescale appearing 
in the problem. We discuss this point below. 



3.1.1. Time step 

The mean rates depend on a as detailed in Paper I (see Fig. 3 of that paper). Rtc{a) is a 
decreasing function of a, and so attains its maximum at a = Omm- AH other rates are either 
constant (ex2), or increasing functions of a, so that their maximum values can be thought 
to occur at a = amax- Accordingly, if we make the following choice for our computational 
time step Atd- 

. f 1 1 1 1 1 1 .^ 

Atd<mm<= ,= ,= ,= , == >, (14) 

K Rtc(^(^min) Rexl(^(^max^ Rex2{S^max) Rdssiflmax) ^couiflmax^ ) 

this will ensure that Atd is smaller than the shortest of the above event-occurrence timescales. 



^However, as is well-known, this time step must also obey the Courant condition (jPress et.al 

19921 ) throughout the range of a under consideration {i.e., Q.QRq-QQRq): 



At, = e^, e < 1. (15) 

J max 

Here, Aa is the step-size in a, and /j^^x is the largest value of /(a) over the range of a 
under co nsideration (see ab ove and paper I). Satisfaction of this condition is essential for the 
stability (jPress et.al Ill992l ) of the solution of Eqn. ([9]). 
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To ensure that both of the above conditions are satisfied, we choose the time step At 
for solving Eqn. iQ to be 

At = min{Atrf, AtJ. (16) 



3.2. Solution of Stochastic Boltzmann Equation 



The Lax-Wendorff scheme ( iPress et.al Ill992l ) used by us for numerical solution of the 
Boltzmann equation in the continuous limit has been introduced in Sec. 2.6 of Paper I. The 
stochastic version of this equation, viz., Eqn. ([9]) can be looked upon as the earlier continuous 
equation with additional stochastic terms, which turns it into a SPDE (see Sec. [T]). We now 
discuss our method of solving this SPDEo. 

It it well-known that ordinary calculus cannot be applied to the handling of stochas- 
tic terms in SPDEs, since these terms are non-differentiable in the ordinary sense, and the 
ordinary definition of an integral does not apply to them. Rather, one has to modify the 
methods of calculus suitably, and redefine appropriate integrals. As summarized in Ap- 
pendix B, one such modified calc ulus is the /to Calculus, which has been used widely for 
solution of SPDEs in recent years (l0ksendal II2004I : iKloeden et.al Ill994j ). The corresponding 
integrals involving the stochastic terms are then called ltd integrals, which have properties 
appropriately different from those of ordinary integrals, as indicated in Appendix B. 



3.2.1. Numerical Method 



In solving an SPDE like Eqn. (Q, one integrates the continu ous terms in the usual way, 
but the stochastic terms must be integrated using Itc3 calculus ( IGains Ill995l ). This means 
that, in advancing the solution at t by a time step dt — which is essentially a Taylor expansion 
of the solution n{a, t) about t — the expansions of the stochastic terms in Eqn. ^ are to be 
performed using the stochastic Taylor expansion (Eqn. (]B7|) ). as discussed in Appendix B. 

A variety of numerical algorithms have been explored by various authors for numer- 
ical solution of SPDEs. The particular algorithm we use is a hybridization of the two- 
step Lax-Wendorff scheme for the continuous terms, as utilized in paper I, and the second- 



^In SPDE literature, the continuous terms are sometimes called drift terms and the stochastic terms 
diffusion terms, but we shall not use this terminology here, since stochastic terms in our problem do not 
always represent diffusion, and furthermore since there is a possibility with such usage of confusion with the 
Fokker-Planck approach, which does represent diffusion in phase space. 
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order stochastic Taylor expans ion according to the Milshtein scheme for the stochastic terms 
(IMilshtein Ill974j : iGains Ill995l ). i.e., Eqn. (|B13p . as explained in Appendix B. In this scheme, 
there is only one stochastic path to be solved for in our case viz., that of n{a, t) (correspond- 
ing to Xfc) and the continuous terms {i.e., the (Tps), the variances in tc, exl, ex2, dss and coll 
rates being as given above. Note that, in each of the two steps in the Lax-Wendorff scheme, 
the expansion (]B13|1 needs to be applied, whereupon we arrive at the following discretization 
schem^ for Eqn. (Q: 



Half step : 



J+1/2 



I (n^i + nf) + R{aj+i/2) - ^(ai+1/2 



„JV , AT 



At 
2 



>N 
+l/2ex'l 



+ 



S„..r, (CLa 



C2 
'-'dss 



N , JV 
2 



2Aa 

Full step 

N+l 

■j 



'i+1 



2Aa 



n 



iV^ 



+ {R{aj) - D{a,)nf) At 



J ex2 J as; 



AT 



+ 



3 ex2' 



S'ex2{^i 



+ m 



dss' 



n 



^dssi^^j 



N 

3 dss 



Aa I '3+1/2 



n 



n 



N 



N+1/2 
J-1/2 



3 coll 

Aa 



n 



N+l/2 
'i+1/2 



2 



^j-1/2 j • 



(17) 



Here, VVj^^ = Sx{0'j)f]^ At, where is the value of a standard normal variate at the A^th 
time step. 



For any particular run, we compute the Wf ^s {Wf^^/^^s) for a particular a 



l'j+l/2) 



over the a and t intervals of integration, and repeat it f or all a,s. The sta ndard normal variate 



rj^'s are generated using the well-known polar method (jPress et.al Ill992l ). All values of Wj ^ 
and y^fj^ii2-^ stored in a two dimensional array {i.e., a Wiener sheet), which serves as 
the input for solving Eqn. (fT7|) . Because of the fluctuations in the coUisional hardening rate 
(as contained in Cacoii)^ impossible that the value of the total hardening rate / might 

occasionally exceed /^axi which would violate the Courant condition, possibly making the 
solution procedure unstable. To avoid this, we have so restricted the variations in VVj^^^^^s 



andW- 



flPress et.alTll992h . 



s th at the amplification factor e = fAt/Aa always lies between zero and unity 



•^It can be shown that the commutation condition (jBlSp is satisfied in this case. 
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4. Results 

We now present the results obtained from our above computations of the cases which 
we studied in Paper 1 in the continuous limit. As before, we study (a) the evolution of the 
distribution function n(a,t), and, (b) the dependence of the computed number of XBs Nxb 
on the Verbunt parameters. We choose exactly the same values of all GC parameters as we 
did in Paper 1, for ease of comparison. 

4.1. Evolution of compact-binary distribution 

In Fig. [H we show a typical evolution of the compact binary population distribution 
n{a,t). The GC parameters were chosen, as in Paper 1, to be p = 6.4 x IO^MqPc"^, = 
0.5 pc and v = 11.6 km sec~^, similar to those of the well-known Galactic cluster 47 Tuc. As 
the figure shows, the surface representing the evolution fluctuates randomly throughout, but 
it does show a clear overall evolution which is of the same nature as that in the continuous 
limit (cf. Fig. 5 of Paper 1). In particular, the population grows with time predominantly at 
shorter radii (a < IORq). As before, we start with a small number of primordial compact 
binaries with various initial distributions, and find that, by t ~ 1 — 1.5 Gyr, the distribution 
"heals" to a form which is independent of the initial choice of distribution. The fluctuations 
differ in detail from run to run, of course, as we choose different seeds for random number 
generation, but the overall nature of the evolution remains the same for all runs. Indeed, the 
results for different runs seem to represent different variations about a mean surface, which 
is very close to that in the continuous limit, as given in Paper 1. We explicitly demonstrate 
this below by displaying temporal and radial slices through the above surface n{a,t) (see 
Paper 1) for different runs, and also displaying their averages over a number of runs, which 
we show to be close to the continuous limit. 

To do this, we first show in Fig. [2] typical time slices, i.e., n{a) at fixed t, (solid lines) 
through the surface in Fig. [H for a single run, overplotting the continuous limit from Paper 
1 for comparison. The distribution with fluctuations does indeed follow the continuous- 
limit distribution generally, the same gross features being visible through fluctuations, in 
particular that n(a) is roughly constant a < 7Rq, and falls off sharply at larger radii. 
The overall nearly-self-similar evolution at large times, described in Paper 1, can also be 
vaguely discerned through the fluctuations. We have discussed possible causes of such self- 
similar evolution in Paper 1. Next, in Fig. [3l we show radial slices corresponding to the 
evolution in Fig. [H representing the behavior of n{t) at a fixed radius a, overplotted with 
the continuous limit from Paper 1. Again, the curves from a single run follow, in a statistical 
sense, the corresponding continuous limits. In particular, it can be seen that the radial slices 
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corresponding to larger values of a tend to saturate by about 6 Gyr, while those for smaller 
values of a do not show such saturation. 

Finally, in Figs. H] and \5[ we show the above temporal and radial slices of the average 
of 12 different runs, overplotted with the the corresponding continuous limits. These figures 
clearly demonstrate how the fluctuations average out over many runs, so that the mean 
result approaches the continuous limit. 



4.2. Number of X-ray binaries 

The total number of GC X-ray binaries Nxb at a particular time was computed from 
Eq. (Q, as in Paper I. We determined Nxb for a representative evolution time of ~ 8 Gyr, and 
studied its dependence on the Verbunt parameters F and 7, so as to relate our computational 
results with the sy stematics of recent observations of X-ray binaries in globular clusters 



( jPooley et.al II2003I ). For this, we computed, as in Paper I, values of Nxb over a rectangular 
grid in F — 7 space, spanning the range 7 = 1 — 10^ and F = 10^ — 10^, which encompasses 
the entire range of Verbunt parameters over which Galactic GCs have been observed (see 
Paper I). Although the GCs actually observed so far lie along a diagonal patch over this grid, 
as explained in Paper I, computational results over the whole grid are useful for clarifying 
the theoretically expected trends. 

As explained in Paper I, at a speciflc grid point (F, 7), the values of p, and Vc are 
evaluated using the deflnitions of Verbunt parameters and the virialization condition (see 
Sec. 3.2 of paper I for a detailed discussion). Also as before, we take representative values of 
primordial stellar binary fraction (kb) and compact-star fraction {kx) to be 10 percent and 
5 percent respectively. 

Fig. [6] shows the resulting Nxb{^,i) surface. As indicated in Paper I, the overall fall- 
off in this surface for 7 > 3 x 10'^ is a signature of the increasing rates of compact-binary 
destruction rates with increasing 7, and the above speciflc value of 7 represents an estimate of 
the threshold above which destruction rates are very important. Further, the trend in Nxb 
with F is simple — Nxb increases with F monotonically, since the dynamical formation 
rate of compact binaries scales with F. What we notice in flg. [6] is that this surface also 
shows random fluctuations due to the stochastic processes, but it generally follows the Nxb 
surface corresponding to the continuous limit, shown overplotted in the same flgure. This 
is similar to what was discussed above for the compact-binary distribution, and the point 
about the mean surface corresponding to the average of many realizations of the stochastic 
processes being very close to the continuous limit also holds here. We also note that the total 
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fluctuations in Nxb increase with increasing value of F. However, as will become evident 
from results discussed below, the relative fluctuations actually decrease with increasing F. 

To further clarify the trends and to make comparisons with the results of the "toy" 
model of B06 and with those of Paper I, we plot the quantity T /Nxb for a fixed value of 
F against 7 in Fig. [71 displaying the curves for several values of F as indicated. As can be 
seen, the fluctuating T /Nxb vs. 7 curves for various values of F follow the same mean trend, 
although the details of the fluctuation are different in different cases. This mean trend is in 
fact very close to the mean "universal" curve corresponding in the continuous limit evolution 
of Paper I, and is overplotted in the figure. Thus, as in the continuous limit case, the basic 
scahng of the toy model, viz., Nxb oc Tg{;~^), where g{;~^) is a "universal" decreasing function 
(representing the increasing binary destruction rate with increasing 7, as explained above), 
does essentially carry over to this detailed model with stochasticity included, suggesting 
a robust feature of the scaling between different clusters which is expected to be further 
confirmed by future observations. 

Another feature of Fig. [7] is that the relative fluctuations in the curves increase with 
decreasing value of F. This is consistent with the intuitive notion that, in all phenomena of 
this nature, the relative fluctuations in Nxb are expected to increase at smaller values of 
NxBi which occur at smaller values of F. More formally, this can be seen as follows. From 
Eqn. (|T3l) . it is clear that, over an interval At, the relative variance in the number of events 
of type X is: 



For the range of F and 7 considered in this work, we found that At was actually close to 
Ate in most cases, so that At ~ 7"^ roughly. Since the formation rates scale as Rx ~ F, we 
have: 



In Sees. 14.11 and 14.21 we saw that the basic trends of the results, as obtained from the 
stochastic Boltzmann equation ([9]), are the same as those obtained from the Boltzmann 
equation in the continuous limit in Paper I. Therefore, as in paper I, the results from the 
stochastic Boltzmann equation are consistent with the observations of XB populations in 
Galactic GCs. Indeed, since fluctuations are present in the dynamical processes under study 
here, we should ideally compare theoretical trends including fluctuations with observational 



rx{a) = (1 -i?x(a)At). 




Therefore, for a flxed 7, rx{a) increases as F (and hence Nxb) decreases. 



4.3. Comparison with observations 
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results, as we do in this paper, where Fig. [6] shows the positions of the observed GCs with 
significant numbers of X-ray sources in the 7 — F — Nxb co-ordinates. The observational 
points do lie near the computed Nxb{.1^'^) surface. In Fig. [71 we compare the T/Nxb — 7 
curves with the positions of the observed points, showing that most points do indeed lie near 
the curves. 

In Fig. [8] we plot the computed contours of constant Nxb in the plane of Verbunt 
parameters, similar to what we did in Fig. 10 of Paper I, but now with the fluctuations 
included. The fluctuations are clearly seen to be larger for smaller values of Nxb, as expected, 
and as mentioned above. Again, the observed numbers generally agree well with the present 
contours which include fluctuations, and these contours do generally follow the continuous- 
limit contours of Paper I, which are shown overplotted. 



Discussions 



We have described in this paper a scheme for introducing stochasticity into the Boltz- 
mann study of compact-binary evolution in globular clusters that we began in Paper I. 
Our scheme involves the use of stochastic calculus (for the first time in this subject, to the 
best of our knowledge), whereas previous studies in the subject have normally used Monte- 



Carlo methods of various descriptions — depending on t 
being studied — for handling s tochasticity (s e e, e.g . 



le particular aspect of the problem 



Hut. McMillan fc Romani I (119921 1: 



Pi Stefano fc Rappaport I (119941 ): iFregeau et.al I (120031 ): iFregeau fc Rasio I (120071 )1 With the 
aid of this scheme, we have demonstrated that the joint action of inherently stochastic and 
continuous processes produces evolutionary trends which necessarily contain fluctuations 
that vary between individual "realizations" of the stochastic processes, as expected. How- 
ever, these trends do generally follow those found in the continuous-limit approximation of 
Paper I, and when trends are averaged over more and more realizations, the mean trend 
comes closer and closer to the continuous-limit one. In this sense, the continuous limit is 
very useful as an indicator of the expected mean trend. On the other hand, the magnitude 
of the fluctuations seen in any given realization, particularly in certain parts of parameter 
space, suggest that one should compare the results of a typical realization to observations, 
in order to get a feel for expected fluctuations in the data from stochastic dynamical pro- 
cesses alone, i.e., apart from those coming from uncertainties in the observational methods 
of obtaining the data. 

Boltzmann approach in its original form appealed to us because of its ability by defini- 
tion to handle weak, frequent, distant encounters and strong, rare, close encounters on the 
same footing. Of course, the approach is of practical use only when probabilities or cross- 
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sections of such encounters are known from detailed studies of individual encounters through 
numerical experiments, as is the case for our current use of this approach. It was generally 
believed that, since Fokker-Planck methods were normally used for handling only the weak, 
frequent, distant encounters above, treating their cumulative effect as a diffusion in phase 
space, this argument would also apply to Monte-Carlo Fokker-Planck methods. However 



in a novel feature inclu ded recently by Fregeau, Rasio and co-authors (IFregeau et.al 112003 



Fregeau fc Rasio 1120071 ) in their Monte-Carlo method, both of the above types of encounters 
are handled in the following way. 

The dynamical evolution of the cluster is treated by a basically Henon-type Monte- 
Carlo method, which describes this evolution as a sequence of equilibrium models, subject 
to regular velocity perturbations which are calculated by the standard Hen on method for 



repres enting the average effect of many weak, frequent, distant encounters (see lFregeau et.al 



(120031 ) and references therein). In addition, the strong, rare, close encounters are by handled 
by (a) keeping track of the (Monte-Carlo-realized) positions of the objects in the cluster, and 
so deciding whether two given objects will undergo a strong, close encounter or not, by a 
rejection method very similar to that described above in Sec. 13.11 and then (b) treating these 
encounters fir st (i) through cross- sections compiled from analytic fits to numerical scattering 
experiments (IFregeau et.al II2003I ). exactly as we have done throughout our approach, and 
then, (ii) in a more detailed approach, through a direct i ntegration of the strong interaction 
at hand using standard two- and three-body integrators ( IFregeau fc Rasio 1120071 ). 



A direct comparison of our results with those of above authors is, for the most part, 
not possible, since we focused primarily on the formation, destruction and hardening of 
a compact binary population in a given GC environment, while Fregeau et. al focused 
primarily on the dynamical evolution of the GC environment in the presence of a given 
primordial binary population. However, there is one feature on which we were able to 
roughly compare our results with those obtained by these and earlier authors. This is the 
problem of hardening of primordial binaries in GCs, pione e ring s tudies which were performed 
through di rect Fokker-Planck integrat i on by iGao et.al I (Il99ll ). and through Monte-Carlo 
method bv lHut. McMillan &: Romani I (119921 ) . and again recently through the above Monte- 
Carlo method by IFregeau et.al I (120031 ). In an early test run of our scheme, we studied this 
problem by "turning off" the binary formation and destruction terms in our scheme, thereby 
studying only the hardening of the primordial binary population through our Boltzmann 
approach. The results we obtained for the progressive hardening of the binary a-distribution 
profile (from an initial profile which was uniform in In a, as in all the above references, and 
in our work) were, indeed, very similar to those given in the above references. 



In a pioneering study, iDi Stefano &: Rappaport I (119921 . Il994j ) explored the tidal-capture 
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formation and subsequent evolution of compact bi naries in GCs, concentrating on recycled 
millisecond pulsars in t he first part of the study ( 



CVs in the second part (jPi Stefano &: Rappaport 



Di Stefano &: Rappaport Ill992l ). and on 
1994h . These authors followed the histo- 



ries of many neutron stars against a given background representing a GC core (parameters 
corresponding to 47 Tuc and u Cen were used as typical examples), employing Monte-Carlo 
methods to generate tidal-capture events in this environment. They followed the subsequent 
orbital evolution of these binaries due to hardening by gravitational radiation and magnetic 
braking, until Roche lobe contact occurred. In those cases where such contact occurred 
through orbit shrinkage before the low-mass companion could reach the giant phase due to 
its nuclear evolution, these authors did not follow further evolution of the binary, while they 
did so when the contact occurred due to the evolutionary expansion of the companion. 

From the above considerations, Di Stefano and Rappaport estimated the expected num- 
ber of recycled pulsars and CVs in GCs like 47 Tuc and uo Cen, and also gave the orbital- 
period distribution of the above binaries at two points, viz., (a) just after tidal capture and 
orbit circularization, and (b) at Roche- lobe contact. However, their orbital-period distri- 
butions cannot be compared directly with those given in this paper (or Paper I) for the 
following reason. In the Monte-Carlo method of these authors, tidal capture occurs at differ- 
ent times for different binaries, as does Roche- lobe contact. Thus, showing the orbital-period 
distribution at any of the above two points means, in effect, that the period-distributions at 
different times are being mixed. By contrast, we have (in this paper and in Paper I) studied 
the evolution of the orbital period-distribution in time, displaying "snapshots" of the whole 
distribution at various times, which we called "time slices" above and in Paper I. In our 
display, for example, at any given time, some binaries are in Roche-lobe contact and some 
are not. Indeed, it seems that the orbit al period-distributions just after tidal capture, as 
given by iDi Stefano fc Rappaport I (119921). sho uld be compared with corresponding N-body 
results given in lPortegies Zwart et.al (1997b), and indeed they appear rather similar. We 
have, of course, pointed out in Paper I, and stress the point here again, that our orbital 
period-distributions are to be regarded at this stage as intermediate steps in our calculation 
— rather than final results to be compared with future data on orbital period-distributions 
of X-ray binaries in GCs — because stellar-evolutionary effects on binary evolution have not 
been included yet in our scheme (also see below). With this inclusion, the aim would be 
to produce the GC -analogue of such orbital period-distributions as have been computed by 
Pfahl et.al. I J^OoJ for LMXBs outside GCs. 



In addition to the above improvement, we listed in Paper I various other improve- 
ments and extensions that are to be implemented in our scheme in future. For example, 
the compact-binary distribution function above can be looked upon as one obtained by in- 
tegrating the full, multivariate distribution function which includes other variables, e.g., the 
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binding energy of the binary in the gravitational potential of the GC — the so-called exter- 
nal binding energy (or, equiyalent ly, the position of the binary within the GC potential well 



( iHut. McMillan fc Romani Ill992l )). over these other variables. It would be most instructive 
to be able to follow the evolution in these additional variables in a more elaborate future 
scheme. 

Encouraged by the veracity of the continuous limit, as presented in this paper, we plan to 
conclude our program of the first stage of exploration of our Boltzmann scheme by studying 
one more problem in the same spirit of demonstration of feasibility as we have followed here 
and in Paper I. This is the question of compact-binary evolution in the environment of an 
evolving GC. Whereas, in keeping with the tradition of numerous previous studies, we have 
treated the GC environment in Paper I and here as a fixed {i.e., unchanging in time) stellar 
background, in reality a GC is believed to undergo considerable evolution following the long, 
quasi-static, "binary-burning" phase, passing through phases of deep core collapse, (possible) 
gravothermal oscillations, and so on. In this concluding study, we propose to demonstrate 
that, at the current level of approximation in our scheme, and in the continuous limit, it is 
possible to follow the evolution of compact-binary populations of GCs through these phases 
of GC evolution, at the expense of only a modest amount of computing time. 

It is a pleasure to thank H. M. Antia, D. Heggie, P. Hut, S. Portegies Zwart, and F. 
Verbunt for stimulating discussions, and the referee for helpful suggestions. 



A. Wiener Processes 

Wiener process is the formal mathematical description of Brownian motion — a classic 
example of a stochastic process, wherein a particle {e.g., pollen grain) on the surface of 
water undergoes random motion due to stochastic bombardment of it by water molecules. A 
standard description of the motion such a particle is given by the following differential form 
due to Langevin: 

dXt = a{t, Xt)dt + a{t, Xt)Ctdt. (Al) 

Here, Xt is one of the components of the particle velocity at time t, and a{t,Xt) is the 
retarding viscous force. The second term on the right-hand side represents the random 
molecular force, represented as a product of an intensity factor a{t,Xt) and a random noise 
factor (t, the latter at each time t being a random number, suitably generated. 

A standard Wiener process W{t) is often defined as a continuous Gaussian process with 
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independent increments, satisfying the following properties: 

W{0) = 0, E{W{t)) = 0, Var(W^(t) - W{s)) = t - s, (A2) 

for all < s < t. Here, E represents the expectation value and 'Var' the variance of the 
indicated"stocLstic variablfi Note that a Wiener proc ess Wt(uj), caii also be thought of 
as a "pure" Brownian motion with a = in Eq. (lAip (jKloeden et.al 1994 ). wherein the 
increments dWt{uj) for any sample path uj represent a Gaussian white noise. 

Eqn. (lAll) can then be rewritten in terms of the symbolic differential (see below) 
dWs^uj) = (s{^^)ds of a Wiener process, and its integral form 

Xt{uj) = Xt,{uj) + j a{s,X,{u))ds+ j a{s,Xs{uj))dWs{u) (A3) 

J U) J to 

represents a path integral over the trajectory of the particle for the sample path Xt{uj), where 
is a particular trajectory of the Brownian particle. 



B. Ito calculus 

The problem with the second term on the right-hand side of Eqn. (lA3p . which represents 
an integral along a Wiener path, is that it is not defined in ordinary calculus, since Wtioj) 
is not differentiable in the ordinary sense. Such an integral along a Wiener path has to be 
redefined suitably to become acceptable mathematically, and the Ito integral is a well-known 
example of this. The classical limit-of-sum definition of an integral does not hold good for 
an Ito integral like 

Xt{uo)= I f{s,uj)dWs{uj), (Bl) 

since the corresponding finite sum will be divergent over a Wiener path, as sample paths of 
a Wiener process do not have bounded variance (see above). Howev er, it can be show n that 



such a sum is mean-square convergent under very general conditions (I0ksendal II2004I ). owing 
to the well-behaved mean-square properties of Wiener processes. Accordingly, Eqn. fIBip is 
defined only in the sense of mean-square convergence, with the result that the integral fIBip 
is a random variable Xtiui) with the following properties: 

E{X,) = 0, E{X^) = f E{f{sf)ds (B2) 



''Strictly speaking, the first equation should be written as 14^(0) 
probability one', since we are dealing with random variables here, 
rigor here, referring the reader to lKloeden et.al J (ll99J) 



= 0, w.p.l, where 'w.p.l' stands for 'with 
. But we shall not go into mathematical 
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Consider now the well-known I to formula for the transformation of a function f{Xt) of a 
stochastic variable Xt ( IGains Ill995l ). For simplicity, first assume that Xf follows a stochastic 
equation of the form 

Xt = Xt^+ [ a{Xt)dt+ [ a{Xt)dWt, (B3) 

Jtn Jtn 



i.e., the same as Eqn. (lA3p . but without explicit time dependence in the continuous and 
stochastic terms. For brevity, we drop the symbol u, representing the sample path, in 
Eqn. flB3|) and from here on. Let us divide the entire time span into time-steps at ti, ^2, • • • 4, • • 
of length hi, h2, . . . hk, . . . with the largest step size hmax- Then Xt at times tk and tk+i are 
related by 



X, 



k+l 



Xk + 



tK+1 



a{Xt)dt 



tK+l 



a{Xt)dWt., 



(B4) 



wher e we write Xk = Xt,. and Xk+i = Xt^^-^ for brevity. The Ito formula states (j0ksendal 



2004f ) that: 



fiXt) = f{Xk) + / Cf{Xs)ds + / f{X,)a{X,)dWs 
Jtk Jtk 

where the operator C is defined by: 



Cf{X,) ^ f'iXMXs) + lnX,)a'iX,). 



(B5) 



(B6) 



For explicitly time-dependent continuous and stochastic terms, the Ito formula can be gen- 
eralized suitably. 

We can use Eqn. f IBSP in Eqn. flB4l) to expand a{Xt) and cr^Xt) around tk'. 



Xk+i = Xk + a{Xk)hk+i + a{Xk)AWk+i 
+ J''^' J' Ca{Xs)dsdt + //;+^ Jl a'{Xs)a{Xs)dWsdt 
+ fi^' d Ca{X,)dsdW, + a'{X.MXs)dWJWt 



(B7) 



'tk •'tk 

Now, if we discard all terms in Eqn. (1B7|) of 0{h°') for a > 1, we obtain 
Xfe+i = Xk + a(X,) + a{Xk)AWk+i + ^a'(X,)a(Xfe) {{AWk+i) 



hk+i) 



:b8) 



which is known as the Milshtein scheme. This is the stochastic analogue of the second-order 
Taylor expansion of ordinary calculus. Th e Milshtein scheme can be shown to be strongly or 
pathwise convergent (IKloeden et.al Ill994l ) to order h, in the sense that the solution converges 
to the actual Brownian path as hmax 0. If we restrict the expansion upto the O^h^^"^) 
terms, i.e., upto the first three terms in the right-hand side of ( ]B7[) . we obtain a slower 
(~ /i^/^) pathwise convergence, which is known as the Euler-Maruyama scheme. 
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For higher dimensions, with Xk. £ "R - ^ and Wt G TV-^ , the second-order stochastic Taylor 
expansion of is given by (see iGains I (119951 ) and references therein): 

Xl^, = Xl + a\X,)h^, + cr]iX,)AW',^, + 5Z S ^^^I^^O^p^I^, k + l) + R, (B9) 

j=l j = l p,q=l 

where ^ ^ 

Ipg{k,k + 1)= [ [ dW'PdW^ (BIO) 

and R contains all terms of 0{h°') for a > 1. If d < p, q {p q) , we obtain upon integration 
by parts: 

I,,{k, k + l) + Igp{k, k + l) = AWj;^,AW',^, = B,g{k, k + 1). (Bll) 
If we further define 

Ap,{k, k+l)= Ipg{k, k+1)- Igp{k, k + 1), (B12) 

then we can, with the aid of Eqns. flB12p and fIBlip . express Ipg in terms of Apg and Bpg. 
Substituting the result in Eqn. ( ]B9p . we finally obtain, 

Xl^, =Xl + a\Xu)h + ^liX,)AWl^, 

p 

j=i p=i 

+ E E 2 b^"^ + a^^0™^-^^'^ + '^ ^^^^^ 

j=l 0<p<q<V ^ ^ 

1 / 9cr* (9cr* \ 



If V 



y f _ = (B14) 

then the terms drop out of Eqn. ( 1B13I) . Equation (1B14I) is called the commutativity 
condition and is usually written as, 

[ap,a,] = 0. (B15) 

When the above commutativity condition is not satisfied, the quantities Apq, known as the 
Levy areas, have to be calculated in order to achieve second-order accuracy. 
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Fig. 1. — A typical example, i.e., one "realization" of the evolution of the binary distribution 
function n{a,t). Globular cluster parameters are chosen to be roughly those of 47 Tuc, as 
explained in text (also see Fig. 5 of Paper I). 
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Fig. 2. — Typical time slices, i.e., n{a) at specified times, for the evolution shown in Fig. [T] 
(solid lines). Overplotted are the same time slices in the continuous limit (dashed lines) from 
Paper I. 
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Fig. 3. — Typical radial slices, i.e., n{t) at fixed values of binary radius for the evolution 
shown in Fig. [H Overplotted are the same radial slices in the continuous limit from Paper I. 
As in Paper I, we show the evolution beyond 8 Gyr by dashed lines to indicate that such long 
evolution times may not be applicable to Galactic GC, but are included here to demonstrate 
the timescales. 
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Fig. 4. — Typical time slices through the average evolutionary surface of 12 different "real- 
izations" of the evolution represented in Fig. [T|, all with the same GC parameters (solid line). 
Overplotted are the corresponding time slices in the continuous limit from paper I (dashed 
line) . 
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Fig. 5. — Typical radial slices of the same average evolutionary surface as in Fig. HI Over- 
plotted are the corresponding radial slices in the continuous limit from paper I. 
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Fig. 6. — Nxb{i, r) surface (solid line). The observed GCs with significant number of XBs 
are shown overplotted. Also shown overplotted is the continuous-limit result from Paper I 
(dashed hne). 
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Fig. 7. — Computed V/Nxb as a function of 7, for values of F as indicated. The continuous- 
limit result for F = 10'' is shown overplotted (thick line). Also shown overplotted are the 
positions of Galctic GCs with significant numbers of X-ray sources, as in Paper I. 
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Fig. 8. — Contours of constant Nxb in the plane of Verbunt parameters. Corresponding con- 
tours in the continuous-hmit case are shown overplotted, using the same hne-styles for easy 
comparison. Positions of GCs with significant numbers of X-ray sources are also overplotted, 
with the corresponding Nxb in parentheses, as in Paper I. 



